function [ CLASSIFICATION ] = mahalanobis( TRAININGDATA, GROUPTRAIN, TESTDATA, modus )
% TRAININGDATA = Matrix mit Trainingsdaten
% GROUPTRAIN = Zuordnung einer Klasse zu den Featurevektoren in TRAINDATA
% TESTDATA = Matrix mit Testdaten
% modus = ist entweder 1,2 or 3
%       modus 1: Berechnen der vollen Kovarianzmatrix fuer jede Klasse
%       modus 2: Annahme, dass die Kovarianzmatrix eine Diagonalmatrix ist
%       modus 3: Annahme, dass alle Klassen die selbe Kovarianzmatrix haben

%% Initialisierung
test_count = size(TESTDATA, 1); % Anzahl der Testdaten
CLASSIFICATION = zeros(test_count, 1); % Initialisierung des Endergebnisses

% Anzahl der vorhandenen Klassen, dabei wird davon ausgegangen, dass die
% Testdaten in GROUPTRAIN mit aufsteigenden Integer repraesentiert werden,
% und alle Elemente einer Klasse in einem zusammenhaengenden Block
% angegeben sind.
% class_count hat so viele Zeilen, wie es Klassen gibt, und enthaelt die
% Anzahl der Elemente jeder Klasse als Eintrag in der zu der Klasse
% korrespondierenden Zeile
class_count = max(max(GROUPTRAIN)); 

% Cell Array, welches die Klassenmittelwerte und Klassenkovarianzen haelt
% wird initialisiert
CLASS_STATS = cell(7, 2); 

% Herausfinden, wieviele objekte jede Klasse hat.
% n(1) sind die Anzahl an Elementen der Klasse 1 usw.
[n, xout] = hist(GROUPTRAIN, 1:class_count);

%% Mittelwerte und Kovarianzmatrizen fuer jede Klasse berechnen (je nach modus)
% Die Matlab version von cov() schaetzt die Kovarianzmatrix mit der 
      %  Annahme, dass der zugrundeliegende Mittelwert auch geschaetzt
      %  wird. Deshalb wird der zusaetzliche Freiheitsgrad in der
      %  Berechnung in den Normalisierungsfaktor einbezogen (1/(N-1). Das
      %  ist genau der richtige Ansatz fuer unsere Berechnung
if modus == 1 % Standard calculation, keine Vereinfachung
    for i = 1:class_count 
        % Der Bereich der aus den Trainingsdaten ausgewaehlt wird, ist
        % immer jener, der alle Featurevektoren einer Klasse beschreibt.
        CLASS_STATS{i,1} = mean( TRAININGDATA( (i-1)*n(i)+1:i*n(i), :) );
        CLASS_STATS{i,2} = cov( TRAININGDATA( (i-1)*n(i)+1:i*n(i), :) );
    end;
elseif modus == 2 % Diagonale Kovarianzmatrizen
    for i = 1:class_count 
        CLASS_STATS{i,1} = mean( TRAININGDATA( (i-1)*n(i)+1:i*n(i), :) );
        CLASS_STATS{i,2} = diag( diag( cov( TRAININGDATA( (i-1)*n(i)+1:i*n(i), :) ) ) );
    end;
elseif modus == 3 % Alle Kovarianzmatrizen sind gleich
    for i = 1:class_count 
        CLASS_STATS{i,1} = mean( TRAININGDATA( (i-1)*n(i)+1:i*n(i), :) );
        CLASS_STATS{i,2} = cov( TRAININGDATA );
    end;
else
    error('Modus muss entweder 1,2 oder 3 sein');
end

%% Klassifikation
% Alle Testobjekte durchlaufen und die Klasse mit der geringsten 
% mahalalobis Distanz fuer das Objekt finden.

for t = 1:test_count
    % Fuer alle Klassen die es gibt
    mahab_class_dist = zeros(7,1);
    for c = 1:class_count
        % Berechnen der Mahalalobis Distanz fuer Klasse c
        p = (TESTDATA(t, :) - CLASS_STATS{c,1});
        q = pinv(CLASS_STATS{c,2}); % Moore-Penrose Invariante kann mit Singularitaeten umgehen.
        u = p';
        mahab_class_dist(c) = sqrt( p * q * u );
    end;
    
    % Auswaehlen der kleinesten Distanz und Zuordnung der passenden Klasse.
    [B, IX] = sort( mahab_class_dist );
    CLASSIFICATION(t) = IX(1); % Der Erste in der Liste ist der mit dem geringsten Abstand.
end;

disp('Mahalanobis Klassifikation abgeschlossen');
